function plotImpulseResponse(IRFy,IRFx,labely,labelx,shockSize,order_app)

fontSizeValue   = 14;
figureSize      =  [0 0 1000 600]; % (0,0) size in x size in y
% Settings for the plotting
defaultFigRow = 3;
defaultFigCol = 3;
normalize     = 0;

% Some dimensions
ne = size(shockSize,1);
nx = size(IRFx,1);
length = size(IRFx,2);

% Name of the shock
tmp = {};
for i=1:size(labelx(1,nx-ne+1:end),2)
    tmp{i} = char(labelx(nx-ne+i));
end
nameOfShocks = tmp(shockSize~=0);

% Removing impulse response function which are zero
[IRFx,labelx] = removeColumn(IRFx,labelx);
[IRFy,labely] = removeColumn(IRFy,labely);

% Merging impulse response functions
IRFall = [IRFy;IRFx];
% Merging labels
labelAll = {};
if ~isempty(labely)
    for i=1:size(labely,2)
        labelAll{i} = char(labely(i));
        % removing potential characters after the symbol '_cu'
        nn = strfind(labelAll{i},char('_cu'));
        if nn > 0
            tmp = labelAll{i};
            labelAll{i} = tmp(1:nn-1);
        end
    end
end
if ~isempty(labelx)
    for i=1:size(labelx,2)
        labelAll{size(labely,2)+i} = char(labelx(i));
        % removing potential characters after the symbol '_cu'
        nn = strfind(labelAll{size(labely,2)+i},char('_cu'));
        if nn > 0
            tmp = labelAll{size(labely,2)+i};
            labelAll{size(labely,2)+i} = tmp(1:nn-1);
        end
    end
end

% Figures and x-axis
totalCharts = size(IRFall,1);
numFig      = ceil(totalCharts/(defaultFigRow*defaultFigCol));
dimX        = 0:size(IRFall,2);

% min and max per impulse response function 
minValues = zeros(1,size(IRFall,1));
maxValues = zeros(1,size(IRFall,1));
for i=1:size(IRFall,1)
    minValues(1,i) = min(min(IRFall(i,:,:)));
    maxValues(1,i) = max(max(IRFall(i,:,:)));
end

idxTotal = 0;
for idxFig=1:numFig
    figure('NumberTitle','off','Name',['Shock to: ' char(nameOfShocks) '; Figure = ' num2str(idxFig)],...
        'Position',figureSize);
    if totalCharts-idxTotal >= 5
        figRow = defaultFigRow;
        figCol = defaultFigCol;
    elseif totalCharts-idxTotal == 4 || totalCharts-idxTotal == 3
        figRow = 2;
        figCol = defaultFigCol;
    else
        figRow = 1;
        figCol = defaultFigCol;
    end
    for i = 1:figRow*figCol
        idxTotal = idxTotal + 1;
        if totalCharts >= idxTotal
            subplot(figRow,figCol,i);
            % Plotting impulse response at first order
            hold on
            if order_app > 1
                % Plotting impulse response at the chosen order
                plot (dimX,[0; reshape(IRFall(idxTotal,:,order_app),length,1)],'-ko');
            end
            plot (dimX,[0; reshape(IRFall(idxTotal,:,1),length,1)],'--k');
            if normalize == 1
                axis([0 length min(minValues) max(maxValues)])
            else
                axis tight
            end
            title(char(labelAll(idxTotal)),'Interpreter','Latex','FontSize',fontSizeValue);
            set(gca,'FontSize',fontSizeValue)
            hold off
        end
        if i == 1
           % Legends
           %if order_app > 1
           %    legend({'Approximation order = 1',['Approximation order = ' num2str(order_app)]},'Orientation','horisontal','Location','NorthEast','EdgeColor',[1 1 1]);
           %else
           %    legend({'Approximation order = 1'},'Orientation','horisontal','Location','NorthEast','EdgeColor',[1 1 1]);
           %end
           %legend({'$\gamma_{\pi} = \hat{\gamma_{\pi}}$','$\gamma_{\pi }=0$'},'Orientation','horisontal','Location','NorthEast','EdgeColor',[1 1 1],...
           %    'Interpreter','Latex','FontSize',fontSizeValue+2);
           legend({'$M^{M,CS}$','$M^{M,CS}$ with $\omega_{d}=0$'},'Orientation','horizontal','Location','NorthEast','EdgeColor',[1 1 1],...
               'Interpreter','Latex','FontSize',fontSizeValue+2);
        end
    end
end
        

% Removing impulse response functions which are zero.
function [IRFout,labelout] = removeColumn(IRFin,labelin)
[n1,n2,n3] = size(IRFin);
labelout = {};
countNoZero = 0;
IRFout     = zeros(countNoZero,n2,n3);
for i=1:n1
    if any(any(abs(squeeze(IRFin(i,:,:)))>1D-10))
        countNoZero = countNoZero + 1;
        IRFout(countNoZero,:,:) = IRFin(i,:,:);
        labelout{countNoZero} = labelin{i};
    end
end

end

    
    
end